#####################################################################################
# Purpose: Create Graphs and Maps
# Code written by: Soeren Henn & Connor Huff
# Last updated: 2023-09-26
#####################################################################################
## This code creates the maps and figures for the paper 
## "The Local Memory of Repression and Who Fights"


rm(list=ls())

## !!!!  set up location of folder here
path <- ""

dir.data <- paste(path, "\\data", sep="")
dir.output <- paste(path, "\\output", sep="")
dir.code <- paste(path, "\\code", sep="")




## libraries
library(lubridate)              # deal with dates
library("stringr")
require("rgdal")                # to ReadOGR
library(rgeos)                  # to gDistance
library(plyr)                   # to join
library(ggplot2)
library(ggsn)       # to add scale bar



#####################################################################################
#################### Figure 1: Population Loss from 1841 to 1851 #################### 
#####################################################################################

########## Panel A: Map
## Load Barony Shapefile
setwd(dir.data)
shape_barony_1841 <- readOGR(dsn = ".", layer = "bar1841")
shape_barony_1841 <- spTransform(shape_barony_1841,
                                 "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +datum=ire65 +units=m +no_defs +ellps=mod_airy +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15")
## Load Covariates data
setwd(dir.data)
data.cov <- read.csv(file="Covariates.csv", stringsAsFactors = F)
names(data.cov)[names(data.cov) == "geog_id"] <- "GEOG_ID"


# merge and plot
shape_barony_1841@data$id = rownames(shape_barony_1841@data)
shape.points = fortify(shape_barony_1841, region="id")
shape.df = join(shape.points, shape_barony_1841@data, by="id")
shape.df2 = join(shape.df, data.cov, by="GEOG_ID")

ggplot(shape.df2, aes(long,lat,group=group, fill = factor(exmortq))) +
  geom_polygon() +
  geom_path(color="black") +
  coord_equal() + 
  north(data=shape.df2, symbol = 16, scale = 0.25, location="topleft") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title=element_blank()) +
  scale_fill_manual(name =element_blank(), 
                    labels = c("Increase", "0-15%", "15-25%", "25-35%", "Over 35%", "No Data"),
                    values = c("grey90", "grey70", "grey50", "grey30", "grey10", "white"))

setwd(dir.output)
ggsave("map_pop_barony_bw.png", width=8, height=8)


########## Panel B: HISTOGRAM OF POPULATION LOSS
setwd(dir.data)
data.merged <- readRDS("Analysis_Data_Barony.rds")

## rural only
data.rural <- data.merged[which(data.merged$rural.1841==1),]

ggplot(data.rural, aes(x=excess.mortality)) +
  geom_histogram(fill="white", color="#009E73", bins=50) +
  theme_bw() +
  coord_cartesian(xlim = c(-.25, .5)) +
  ylab("Count") +
  xlab("Population Loss 1841-51 (in %)")
setwd(dir.output)
ggsave("fig_pop_loss_hist.png", width=6, height=6)


####################################################################################
################# Figure 2: Conflict Participation Rates by Barony ################# 
####################################################################################


setwd(dir.analysis)
data.merged <- readRDS("Analysis_Data_Barony.rds")

## create quartiles
# militia
summary(data.merged$militia.pc.all01)
data.merged$militia.q <- NA
data.merged$militia.q[which(data.merged$militia.pc.all01<0.005)] <- 1
data.merged$militia.q[which(data.merged$militia.pc.all01>0.005 & data.merged$militia.pc.all01<0.013)] <- 2
data.merged$militia.q[which(data.merged$militia.pc.all01>0.013 & data.merged$militia.pc.all01<0.033)] <- 3
data.merged$militia.q[which(data.merged$militia.pc.all01>0.033)] <- 4
data.merged$militia.q[which(data.merged$rural==0)] <- 0
table(data.merged$militia.q)
# WW1 enlistment
summary(data.merged$enlist.pc.all)
data.merged$enlist.q <- NA
data.merged$enlist.q[which(data.merged$enlist.pc.all<0.006)] <- 1
data.merged$enlist.q[which(data.merged$enlist.pc.all>0.006 & data.merged$enlist.pc.all<0.01)] <- 2
data.merged$enlist.q[which(data.merged$enlist.pc.all>0.01 & data.merged$enlist.pc.all<0.018)] <- 3
data.merged$enlist.q[which(data.merged$enlist.pc.all>0.018)] <- 4
data.merged$enlist.q[which(data.merged$rural==0)] <- 0
table(data.merged$enlist.q)
data.merged2 <- subset(data.merged, select=c(GEOG_ID, militia.q, enlist.q))
# WW1 casualties
summary(data.merged$deaths.pc.all)
data.merged$death.q <- NA
data.merged$death.q[which(data.merged$deaths.pc.all<0.0045)] <- 1
data.merged$death.q[which(data.merged$deaths.pc.all>0.0045 & data.merged$deaths.pc.all<0.0073)] <- 2
data.merged$death.q[which(data.merged$deaths.pc.all>0.0073 & data.merged$deaths.pc.all<0.013)] <- 3
data.merged$death.q[which(data.merged$deaths.pc.all>0.013)] <- 4
data.merged$death.q[which(data.merged$rural==0)] <- 0
table(data.merged$death.q)
data.merged2 <- subset(data.merged, select=c(GEOG_ID, militia.q, enlist.q, death.q))
# IRA
summary(data.merged$IRA.pc.all)
data.merged$IRA.q <- NA
data.merged$IRA.q[which(data.merged$IRA.pc.all<0.001272)] <- 1
data.merged$IRA.q[which(data.merged$IRA.pc.all>0.001272 & data.merged$IRA.pc.all<0.002398)] <- 2
data.merged$IRA.q[which(data.merged$IRA.pc.all>0.002398 & data.merged$IRA.pc.all<0.003961)] <- 3
data.merged$IRA.q[which(data.merged$IRA.pc.all>0.003961)] <- 4
data.merged$IRA.q[which(data.merged$rural==0)] <- 0
table(data.merged$IRA.q)
data.merged2 <- subset(data.merged, select=c(GEOG_ID, militia.q, enlist.q, death.q, IRA.q))



## read shapefile
setwd(dir.data)
shape_barony_1841 <- readOGR(dsn = ".", layer = "bar1841")
shape_barony_1841 <- spTransform(shape_barony_1841,
                                 "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +datum=ire65 +units=m +no_defs +ellps=mod_airy +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15")


shape_barony_1841@data$id = rownames(shape_barony_1841@data)
shape.points = fortify(shape_barony_1841, region="id")
shape.df = join(shape.points, shape_barony_1841@data, by="id")
shape.df2 = join(shape.df, data.merged2, by="GEOG_ID")

ggplot(shape.df2, aes(long,lat,group=group, fill = factor(militia.q))) +
  geom_polygon() +
  geom_path(color="black") +
  coord_equal() + 
  north(data=shape.df2, symbol = 16, scale = 0.25, location="topleft") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title=element_blank()) +
  scale_fill_manual(name =element_blank(), 
                    labels = c("Urban", "0-0.5%", "0.5-1.3%", "1.3-3.3%", "Over 3.3%", "No Data"),
                    values = c("grey90", "grey70", "grey50", "grey30", "grey10", "white"))



setwd(dir.output)
ggsave("map_conflict_militia_bw.png", width=8, height=8)

ggplot(shape.df2, aes(long,lat,group=group, fill = factor(enlist.q))) +
  geom_polygon() +
  geom_path(color="black") +
  coord_equal() + 
  north(data=shape.df2, symbol = 16, scale = 0.25, location="topleft") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title=element_blank()) +
  scale_fill_manual(name =element_blank(), 
                    labels = c("Urban", "0-0.6%", "0.6-1%", "1-1.8%", "Over 1.8%", "No Data"),
                    values = c("grey90", "grey70", "grey50", "grey30", "grey10", "white"))

setwd(dir.output)
ggsave("map_conflict_enlist_bw.png", width=8, height=8)


ggplot(shape.df2, aes(long,lat,group=group, fill = factor(death.q))) +
  geom_polygon() +
  geom_path(color="black") +
  coord_equal() + 
  north(data=shape.df2, symbol = 16, scale = 0.25, location="topleft") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title=element_blank()) +
  scale_fill_manual(name =element_blank(), 
                    labels = c("Urban", "0-0.45%", "0.45-0.73%", "0.73-1.3%", "Over 1.3%", "No Data"),
                    values = c("grey90", "grey70", "grey50", "grey30", "grey10", "white"))

setwd(dir.output)
ggsave("map_conflict_death_bw.png", width=8, height=8)



ggplot(shape.df2, aes(long,lat,group=group, fill = factor(IRA.q))) +
  geom_polygon() +
  geom_path(color="black") +
  coord_equal() +  
  north(data=shape.df2, symbol = 16, scale = 0.25, location="topleft") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title=element_blank()) +
  scale_fill_manual(name =element_blank(), 
                    labels = c("Urban", "0-0.13%", "0.13-0.24%", "0.24-0.4%", "Over 0.4%", "No Data"),
                    values = c("grey90", "grey70", "grey50", "grey30", "grey10", "white"))

setwd(dir.output)
ggsave("map_conflict_IRA_bw.png", width=8, height=8)




####################################################################################
################################ Figure 3: Raw Data ################################
####################################################################################
setwd(dir.analysis)
data.merged <- readRDS("Analysis_Data_Barony.rds")

## rural only
data.rural <- data.merged[which(data.merged$rural.1841==1),]
summary(data.rural$enlist.pc.all)

## population loss only
data.rural <- data.rural[which(data.rural$excess.mortality>=0),]

## Enlistment WW1
ggplot(data.rural, aes(x=excess.mortality,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Enlistment WW1 p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_enlistment_loss.png", width=6, height=6)

## Casualties WW1
ggplot(data.rural, aes(x=excess.mortality,y=deaths.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Casualties WW1 p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_death_loss.png", width=6, height=6)

ggplot(data.rural, aes(x=excess.mortality,y=militia.pc.all01)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Militia members p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_Militia_loss.png", width=6, height=6)

ggplot(data.rural, aes(x=excess.mortality,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Irish Rebels p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_IRA_loss.png", width=6, height=6)


##############################################################################
#### Figure 4: Correlation Sinn Fein Vote Share and Conflict Participation ###
##############################################################################

setwd(dir.data)
data.rural <- readRDS("Analysis_Data_Const.rds")

## Enlistment WW1
ggplot(data.rural, aes(x=share.SF.max,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Enlistment WW1 p/c") +
  xlab("Sinn Fein Vote Share 1918")
setwd(dir.output)
ggsave("fig_raw_enlistment_SF.png", width=6, height=6)

ggplot(data.rural, aes(x=share.SF.max,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Irish Rebels p/c") +
  xlab("Sinn Fein Vote Share 1918")
setwd(dir.output)
ggsave("fig_raw_IRA_SF.png", width=6, height=6)



##############################################################################
## Figure 5: Correlation Between 1811 Economic Indicators and Conflict Participation ##
##############################################################################

setwd(dir.data)
data.merged <- readRDS("Analysis_Data_Barony.rds")

## rural only
data.rural <- data.merged[which(data.merged$rural.1841==1),]
summary(data.rural$enlist.pc.all)

## only keep baronies with population loss
data.rural <- data.rural[which(data.rural$excess.mortality>=0),]

########## Panel A: Enlistment WW1 & labourer 1911
ggplot(data.rural, aes(x=labourer2.perc.1911,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Enlistment WW1 p/c") +
  xlab("Percentage Labourer 1911")
setwd(dir.output)
ggsave("fig_raw_enlistment_lab11.png", width=6, height=6)

########## Panel B: Enlistment WW1 & literacy 1911
ggplot(data.rural, aes(x=read.write.1911,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Enlistment WW1 p/c") +
  xlab("Literacy 1911")
setwd(dir.output)
ggsave("fig_raw_enlistment_lit11.png", width=6, height=6)


########## Panel C: IRA & labourer 1911
ggplot(data.rural, aes(x=labourer2.perc.1911,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Irish Rebels p/c") +
  xlab("Percentage Labourer 1911")
setwd(dir.output)
ggsave("fig_raw_IRA_lab11.png", width=6, height=6)

########## Panel D: IRA & literacy 1911
ggplot(data.rural, aes(x=read.write.1911,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Irish Rebels p/c") +
  xlab("Literacy 1911")
setwd(dir.output)
ggsave("fig_raw_IRA_lit11.png", width=6, height=6)




##############################################################################
## Figure A1: Correlation Between 1911 Literacy and Conflict Participation by Religion ##
##############################################################################

setwd(dir.data)
data.merged <- readRDS("Analysis_Data_Barony.rds")

## rural only
data.rural <- data.merged[which(data.merged$rural.1841==1),]
summary(data.rural$enlist.pc.all)

## only keep baronies with population loss
data.rural <- data.rural[which(data.rural$excess.mortality>=0),]

########## Panel A: Enlistment WW1 & Catholic literacy
ggplot(data.rural, aes(x=literacy.catholic.1911,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Enlistment WW1 p/c") +
  xlab("Catholic Literacy 1911")
setwd(dir.output)
ggsave("fig_raw_enlistment_lit11c.png", width=6, height=6)

########## Panel B: IRA & Catholic literacy
ggplot(data.rural, aes(x=literacy.catholic.1911,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Irish Rebels p/c") +
  xlab("Catholic Literacy 1911")
setwd(dir.output)
ggsave("fig_raw_IRA_lit11c.png", width=6, height=6)

########## Panel C: Enlistment WW1 & Non-Catholic literacy
ggplot(data.rural, aes(x=literacy.noncatholic.1911,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Enlistment WW1 p/c") +
  xlab("Non-Catholic Literacy 1911")
setwd(dir.output)
ggsave("fig_raw_enlistment_lit11nc.png", width=6, height=6)

########## Panel D: IRA & Non-Catholic literacy
ggplot(data.rural, aes(x=literacy.noncatholic.1911,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  ylab("Irish Rebels p/c") +
  xlab("Non-Catholic Literacy 1911")
setwd(dir.output)
ggsave("fig_raw_IRA_lit11nc.png", width=6, height=6)

##############################################################################
################## Figure A2: Raw Data with Population Gain ################## 
##############################################################################

setwd(dir.data)
data.merged <- readRDS("Analysis_Data_Barony.rds")

## rural only
data.rural <- data.merged[which(data.merged$rural.1841==1),]
summary(data.rural$enlist.pc.all)

## fix one outlier
data.rural$excess.mortality[which(data.rural$excess.mortality< -0.18)] <- -0.18

########## Panel A: Militia Members
ggplot(data.rural, aes(x=excess.mortality,y=militia.pc.all01)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  #coord_cartesian(ylim = c(-1.5, 2)) +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Militia members p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_Militia.png", width=6, height=6)

########## Panel B: Enlistment WW1
ggplot(data.rural, aes(x=excess.mortality,y=enlist.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Enlistment WW1 p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_enlistment.png", width=6, height=6)

########## Panel C: Casualties WW1
ggplot(data.rural, aes(x=excess.mortality,y=deaths.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Casualties WW1 p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_death.png", width=6, height=6)

########## Panel D: Irish Rebels
ggplot(data.rural, aes(x=excess.mortality,y=IRA.pc.all)) +
  geom_point(data=data.rural,color="grey",size=2,alpha = .4) +
  stat_summary_bin(fun.y='mean', bins=25,
                   color="#009E73", size=2, geom='point', data=data.rural) +
  geom_smooth(method = 'lm', color='black', data=data.rural) +
  theme_bw() +
  #coord_cartesian(ylim = c(-1.5, 2)) +
  geom_vline(xintercept = 0, linetype="dotted", color = "black", size=1) +
  ylab("Irish Rebels p/c") +
  xlab("Population Loss 1841-51")
setwd(dir.output)
ggsave("fig_raw_IRA.png", width=6, height=6)






#####################################################################################
######################## Figure A3: Sinn Fein Vote Share 1918 ####################### 
#####################################################################################

########## Election Results
setwd(dir.data)
data.analysis <- readRDS("Analysis_Data_Const.rds")

summary(data.analysis$share.SF.max)
data.analysis$SF.q <- NA
data.analysis$SF.q[which(data.analysis$share.SF.max<0.50)] <- 1
data.analysis$SF.q[which(data.analysis$share.SF.max>0.50 & data.analysis$share.SF.max<0.70)] <- 2
data.analysis$SF.q[which(data.analysis$share.SF.max>0.70 & data.analysis$share.SF.max<0.869)] <- 3
data.analysis$SF.q[which(data.analysis$share.SF.max>0.869)] <- 4
data.analysis$SF.q[which(data.analysis$rural==0)] <- 0
table(data.analysis$SF.q)
data.analysis$CON_NAME <- data.analysis$Con_name

########## Constituency Boundaries
setwd(dir.data)
shape.election <- readOGR(dsn = ".", layer = "constituencies2")
shape.election <- spTransform(shape.election,
                                 "+proj=tmerc +lat_0=53.5 +lon_0=-8 +k=1.000035 +x_0=200000 +y_0=250000 +datum=ire65 +units=m +no_defs +ellps=mod_airy +towgs84=482.530,-130.596,564.557,-1.042,-0.214,-0.631,8.15")


shape.election@data$id = rownames(shape.election@data)
shape.points = fortify(shape.election, region="id")
shape.df = join(shape.points, shape.election@data, by="id")
shape.df2 = join(shape.df, data.analysis, by="CON_NAME")

########## Panel A: Map
ggplot(shape.df2, aes(long,lat,group=group, fill = factor(SF.q))) +
  geom_polygon() +
  geom_path(color="black") +
  coord_equal() + 
  north(data=shape.df2, symbol = 16, scale = 0.25, location="topleft") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(),
        axis.line = element_blank(), axis.text = element_blank(), axis.ticks = element_blank(), axis.title=element_blank()) +
  scale_fill_manual(name =element_blank(), 
                    labels = c("0-50%", "50-70%", "Over 70%", "Unopposed", "No Data"),
                    values = c("grey80", "grey60", "grey40", "grey20", "white"))
setwd(dir.output)
ggsave("map_Sinn_Fein_bw.png", width=8, height=8)

########## Panel B: Histogram
ggplot(data.analysis, aes(x=share.SF.max)) +
  geom_histogram(fill="white", color="#009E73", bins=50) +
  theme_bw() +
  #coord_cartesian(xlim = c(-.25, .5)) +
  ylab("Count") +
  xlab("1918 Sinn Fein Vote Share")
setwd(dir.output)
ggsave("fig_Sinn_Fein_hist.png", width=6, height=6)


##################################################################
################# Figure A5: Sensitivity Analysis ################ 
##################################################################

setwd(dir.data)

data.merged <- readRDS("Analysis_Data_Barony.rds")

## rural only
data.rural <- data.merged[which(data.merged$rural.1841==1),]
summary(data.rural$enlist.pc.all)

## only keep baronies with population loss
data.rural <- data.rural[which(data.rural$excess.mortality>=0),]

reg.enlist <- lm(enlist.pc.all ~ excess.mortality + geo.distbattles + perc.cath.1841 + perc.agri.1841 + log.area41 + pop.den.1841 + literacy.1841 + log.pop41 + housing.fourth.1841 + potato.suit.bar + july.temperature + july.rainfall +  geo.rugged + geo.distcoast + geo.distBelfast + geo.distDublin + factor(county_id), data=data.rural, na.action=na.omit)
summary(reg.enlist)
enlist.sensitivity <- sensemakr(model = reg.enlist, 
                                treatment = "excess.mortality",
                                benchmark_covariates = "perc.cath.1841",
                                kd = c(5,15,25),
                                #ky = 1:3, 
                                q = 1,
                                alpha = 0.05, 
                                reduce = TRUE)
summary(enlist.sensitivity)
plot(enlist.sensitivity)
setwd(dir.output)
pdf("fig_sensitivity_enlist.pdf")
plot(enlist.sensitivity)
dev.off()

## casualties
reg.death <- lm(deaths.pc.all ~ excess.mortality + geo.distbattles + perc.cath.1841 + perc.agri.1841 + log.area41 + pop.den.1841 + literacy.1841 + log.pop41 + housing.fourth.1841 + potato.suit.bar + july.temperature + july.rainfall +  geo.rugged + geo.distcoast + geo.distBelfast + geo.distDublin + factor(county_id), data=data.rural, na.action=na.omit)

death.sensitivity <- sensemakr(model = reg.death, 
                               treatment = "excess.mortality",
                               benchmark_covariates = "perc.cath.1841",
                               kd = c(5,15,25),
                               #ky = 1:3, 
                               q = 1,
                               alpha = 0.05, 
                               reduce = TRUE)
summary(death.sensitivity)
setwd(dir.output)
pdf("fig_sensitivity_death.pdf")
plot(death.sensitivity)
dev.off()

## Militia
reg.militia <- lm(militia.pc.post01 ~ excess.mortality + geo.distbattles + perc.cath.1841 + perc.agri.1841 + log.area41 + pop.den.1841 + literacy.1841 + log.pop41 + housing.fourth.1841 + potato.suit.bar + july.temperature + july.rainfall +  geo.rugged + geo.distcoast + geo.distBelfast + geo.distDublin + factor(county_id), data=data.rural, na.action=na.omit)

militia.sensitivity <- sensemakr(model = reg.militia, 
                                 treatment = "excess.mortality",
                                 benchmark_covariates = "perc.cath.1841",
                                 kd = c(5,15,25,50),
                                 #ky = 1:3, 
                                 q = 1,
                                 alpha = 0.05, 
                                 reduce = TRUE)
summary(militia.sensitivity)
setwd(dir.output)
pdf("fig_sensitivity_militia.pdf")
plot(militia.sensitivity)
dev.off()

## IRA

reg.IRA <- lm(IRA.pc.all ~ excess.mortality + geo.distbattles + perc.cath.1841 + perc.agri.1841 + log.area41 + pop.den.1841 + literacy.1841 + log.pop41 + housing.fourth.1841 + potato.suit.bar + july.temperature + july.rainfall +  geo.rugged + geo.distcoast + geo.distBelfast + geo.distDublin + factor(county_id), data=data.rural, na.action=na.omit)

IRA.sensitivity <- sensemakr(model = reg.IRA, 
                             treatment = "excess.mortality",
                             benchmark_covariates = "perc.cath.1841",
                             kd = c(5,10,15),
                             #ky = 1:3, 
                             q = 1,
                             alpha = 0.05, 
                             reduce = TRUE)
summary(IRA.sensitivity)
plot(IRA.sensitivity)
setwd(dir.output)
pdf("fig_sensitivity_IRA.pdf")
plot(IRA.sensitivity)
dev.off()


